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KHOVANSKII-ROLLE CONTINUATION FOR REAL SOLUTIONS 

DANIEL J. BATES AND FRANK SOTTILE 

Abstract. We present a new continuation algorithm to find all nondegenerate real solu- 
tions to a system of polynomial equations. Unlike homotopy methods, it is not based on 
a deformation of the system; instead, it traces real curves connecting the solutions of one 
system of equations to those of another, eventually leading to the desired real solutions. 
It also differs from homotopy methods in that it follows only real paths and computes no 
complex solutions of the original equations. The number of curves traced is bounded by 
the fewnomial bound for real solutions, and the method takes advantage of any slack in 
that bound. 



Introduction 

Numerical continuation gives efficient algorithms for finding all complex solutions to a 
system of polynomial equations jlSj. Current implementations |H |15], [19] routinely and 
reliably solve systems with thousands of solutions having dozens of variables. Often it is 
the real solutions or the solutions with positive coordinates that are sought. In these cases, 
numerical continuation first finds all complex solutions, which are then sifted to obtain 
the desired real ones, i.e., those with sufficiently small imaginary part. This is inefficient, 
particularly for fewnomial systems (systems of polynomials with few monomials [13]) which 
possess bounds on their numbers of positive [7] and of real [3] solutions that can be far 
smaller than their number of complex solutions. More fundamentally, homotopy methods 
do not exploit the real algebraic nature of the problem. 

We present a numerical continuation algorithm to find all nondegenerate real solutions 
(or all positive solutions) to a system of polynomials. This Khovanskii-RoUe continuation 
algorithm is modeled on the derivation of fewnomial bounds [3l [7] for the numbers of real 
solutions. It is efficient in that the number of paths followed depends upon the correspond- 
ing fewnomial bounds and not on the number of complex solutions. This algorithm also 
naturally takes advantage of any slack in the fewnomial bound. All path tracing is along 
real curves and except for a precomputation, the algorithm is completely real. 

To the best of our knowledge, this is the first continuation algorithm that finds only 
real solutions to a system of equations. There are other numerical algorithms that find 
only real solutions. Exclusion [18], §6.1] recursively subdivides a bounded domain exclud- 
ing subdomains that cannot harbor solutions. This is an active area with research into 
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better exclusion tests [TOl [16] and complexity [T]. Lasserre, Laurent, and Rostalski [T^ 
recently proposed another method based on sums-of-square and semidefinite programming. 
Cylindrical algebraic decomposition [9j| is a symbolic method that computes of a detailed 
stratification of the ambient space in addition to the solutions. Grobner bases may be used 
to compute a rational univariate representation [17] from which the solutions are found by 
substituting the roots of a univariate polynomial into the representation. 

Many existing methods have drawbacks. Exclusion suffers from the curse of dimension- 
ality: the growth in the number of cells may be exponential in the number of variables. 
Cylindrical algebraic decomposition is infeasible in dimensions above 3 due to the com- 
plexity of the data structures involved. Grobner basis algorithms essentially compute all 
complex solutions and much more in addition to the real solutions. 

Our algorithm is based on Khovanskii's generalization of RoUe's Theorem |T3] and the 
proof of the fewnomial bound [3l [7] , and we are indebted to Bernd Sturmfels who observed 
that this proof leads to a numerical algorithm to compute real solutions to systems of 
polynomials. This algorithm does not directly solve the given polynomial system but rather 
an equivalent (Gale dual |8j) system consisting of master functions in the complement of 
an arrangement of hyperplanes. If the original system involves n polynomials in n variables 
having a total of l+n+1 monomials, then the Gale dual system consists of / master functions 
in the complement of an arrangement of l+n+1 hyperplanes in R'. Gale duality [8] asserts 
that the two systems are equivalent and produces an algebraic bijection between complex 
solutions, which restricts to a bijection between their real solutions, and further restricts to 
a bijection between their (appropriately defined) positive solutions. 

Here, we describe the Khovanskii-Rolle algorithm and a proof-of-concept implementation. 
This implementation uses the system Bertini j|| for precomputation, but it is otherwise 
implemented in Maple. It is also restricted to the case of finding positive solutions when 
/ = 2. Detailed examples are on our webpage [6], which also contains our software and 
a brief guide. We plan significant further work on this algorithm, including source code 
and compiled binaries of an implementation that places no restrictions on / and links to 
numerical homotopy continuation software [U |15l |T9] for the precomputation. We also plan 
theoretical work on the technical challenges posed by this algorithm. 

As the algorithm does not solve the original system, but rather an equivalent system in a 
different number of variables, we first explain that transformation, together with an example 
and an important reduction in Section [H We also discuss curve tracing and homotopy 
continuation in Section [H In Section we explain our algorithm, give a proof of correctness, 
and examine its complexity. We illustrate the Khovanskii-Rolle algorithm on an example 
in Section [3] and describe the performance of our implementation on this example and on 
one that is a bit more extreme. Surprisingly, our simple maple implementation for finding 
positive solutions outperforms both Bertini |4j and PHCpack ^J] . In Section HI we describe 
our implementation, and we end the paper with a discussion of future work. Please note 
that we intend "nondegenerate real solutions" when we say "real solutions" throughout, 
unless otherwise noted. 
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1. Background 

We first explain liow Gale duality transforms a system of polynomial equations into an 
equivalent system of rational master functions. This reduces the problem of finding all 
positive solutions to the polynomial equations to finding solutions to the master functions 
in a polyhedron, which we may assume is bounded. Then we discuss general curve-tracing 
algorithms and finally outline the method of homotopy continuation to highlight how it 
differs from Khovanskii-RoUe continuation. 



1.1. Gale Duality. The Khovanskii-Rolle algorithm does not in fact directly solve polyno- 
mial systems, but rather solves systems of master functions defined in the complement of 
an arrangement of hyperplanes. Solutions of the system of master functions correspond to 
solutions of the original system via Gale duality [S]. We begin with an example. 



Example 1.1. The system of Laurent polynomials (written in diagonal form) 
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we have the following identity on the monomials 
(1.2) 



cd)\be^)-\cd-^e-^)\bc-^e 
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A Gale system dual to (11. ip is obtained from the identity (II. 2p by first substituting x for 
the term jbe^ and y for a~^b~^e in (II. ip to obtain 



cd 
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3-2x + y , 



cd^^e^^ 
bc~'^e 
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Then, we substitute these hnear polynomials for the monomials in (11.21) to obtain the system 
of rational master functions 

(1.3) ^ \^ % 

y(2x+2y-ir(4x)-3(i±f^) (e^j (^8:^) (3_2x+y) = 1. 

If we solve these for 0, they become f = g = 0, where 

/ := (2x+2y-l)(l+x-y)(8-3x-2?/)2 - 8yx'^(6-x-3yy(3-2x+y) , 

(1.4) 

g := y{2x+2y-l)%l+x~y)%8-3x-2yy{3-2x+y) - 327Q8x^{6-x-3yf . 

Figure [T] shows the curves defined by / and g and the lines given by the linear factors in 
/ and g. The curve / = has the three branches indicated and the other arcs belong to 
^ = 0. 




3-2x+y = 



Figure 1. Curves and lines 

It is clear that the solutions to f = g = in the complement of the lines are consequences 
of solutions to (II. ip . More, however, is true. The two systems define isomorphic schemes 
as complex or as real varieties, with the positive solutions to (11.11) corresponding to the 
solutions of / = (7 = lying in the central heptagon. Indeed, the polynomial system (II. ip 
has 102 solutions in (C^)^. This may be verified with either symbolic software such as 
Singular [11], or with numerical solvers. Ten of these solutions are real, with six positive. 
We computed them using both Bertini [1] and PHCpack [19] (which produce the same 
solutions) and display them (in coordinates {a,b,c,d,e)) to 3 significant digits. 

(-8.92, -1.97, -0.690, 3.98, -1.28) , (-0.0311, 8.52, -1.26, -14.1, -1.39) , 

(-0.945, 3.41, 1.40, 0.762, 1.30) , (-5.21, 4.57, 2.48, 1.20, 1.34) , 

(1.5) (2.19, 0.652, 0.540, 2.27, 1.24) , (2.45, 0.815, 0.576, 1.38, 1.20) , 

(3.13, 1.64, 0.874, 0.962, 1.28) , (2.00, 0.713, 1.17, 3.28, 2.20) , 

(1.61, 1.03, 2.37, 1.98, 2.35) , (0.752, 3.10, 2.36, 1.55, 1.48) . 
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For the system / = (7 = of master functions (11.41) in the complement of the hnes shown 
in Figure [H symbohc and numerical software finds 102 solutions with ten real solutions and 
six lying in the heptagon of Figure [H We give the real solutions to 3 significant digits 

(-0.800,-0.0726), (4.13,5.26), (1.44,-0.402), (2.04,-0.0561), 
(0.249, 0.864), (0.296, 0.602), (0.670,0.250), (0.866,1.54), (1.43,1.42), (1.70,0.634). 

These correspond (in order) to the solutions in the list (11. 5p under the map (a, b, c, d, e) 1— > 
(6e^/4, ea~^6~^). Gale duality generalizes this equivalence. 

Let A = {0, C be a collection of 1+l+n integer vectors, which we 
consider to be exponents of (Laurent) monomials: If a = (ai, . . . , a„), then 
A polynomial with support .4. is a linear combination of monomials from A, 

f ■ ~ ^ ^ Cq, X , 

where Cq G M - our poljTiomials and functions are real. Since we allow negative exponents, 
/ is defined on the non-zero real numbers, (M^)*^. We consider systems of polynomials 

(1.6) fi{Xi,...,Xn) = /2(Xi, . . . ,X„) = ■■■ = fn{Xi,...,Xn) = 0, 

in which each polynomial has the same support A. We also assume that the system is 
general, by which we mean that it defines a finite set of non-degenerate points in (M^)", in 
that the differentials of the polynomials fi are linearly independent at each solution. 

This implies that the polynomials /, are linearly independent over M, in particular that 
the coefficient matrix of n of the monomials, say x"^, . . . ,a;°", is invertible. We may then 
solve (II. 6p for these monomials to obtain a diagonal system 

(1.7) x"' = pi(a;""+\...,x°'+"), for z = 1, . . . , n , 

where the pt are degree 1 polynomials in their arguments. This has the same form as (11. II) . 

Let B := {(3i, . . . , Pi} C Z'^" be a basis for the free Abelian group A'^ of integer linear 
relations among the vectors in A, 

A^ := {(5i, . . . , k+n) G Z'+" I biai + ■■■ + hi+nai+n = 0} . 

If we write (3j = {jSj^i \ i = 1, . . . , l+n), then the monomials x"' satisfy 

1 = (a;°i)^^.i...(x"'+")^^-'+", for J = 1,...,/. 

Using the diagonal system (11.71) to substitute Pi = pi{x°'"+^, . . . ,x"'+") for we obtain 

1 = -1 . . .pl^j^" . . . . ^ for J = 1, . . . , / . 

If we now set i/i := x°'"+' and write pn+i{y) = for i = 1, . . . , / , then we get 

(1.8) 1 = Piiyf-'---Pi+niyf^-'^-, forj = l,...,/. 

These polynomials Piiy) define an arrangement Ti, of hyperplanes in either W or C'. 
Because the exponents (3ij may be negative, the system (11.81) only makes sense in the 
complement M-j-i of this arrangement. These rational functions, Pi{y)^^'^ ■ ■ -pi+niy)^^'^^", are 
called master functions, and they arise in the theory of hyperplane arrangements. We give 
the main result of [8]. 
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Theorem 1.2. If the integer linear span lA oj the exponents A is equal to Z", then the as- 
sociation ^ Hi defines a bijection between solutions to the system of polynomials (11. 6p 
in (C^)" and solutions to the system of master functions (11. 8p in M-ji^. 

This bijection restricts to a bijection between real solutions to both systems and further to 
a bijection between the positive solutions to (11.61) and the solutions to the system of master 
functions (II. 8p which lie in the positive chamber A := {y & M! \ Pi{y) > 0, i = 1, . . . ,1 + n} . 

The bijections of Theorem 11.21 respect the muhiphcities of the solutions. Because of this 
equivalence between the two systems, and because the exponents A and B come from Gale 
dual vector configurations, we call the system (II. 8p of master functions Gale dual to the 
original polynomial system. 

Theorem 11.21 may be strengthened as follows. There still is a bijection between real 
solutions if we only require that ZA have odd index in Z"' and that the span of . . . , 
have odd index in A'^. The bijection between positive solutions in M" and solutions in the 
positive chamber A in M' remains valid if we allow A and Pi to be any real vectors such that 
A spans M" and . . . , Pi} is a basis of the real vector space A-^. The Khovanskii-Rolle 
continuation algorithm is naturally formulated in terms of real exponents Pi. 

1.1.1. Master functions on bounded poly hedra. Before we discuss the numerical algorithms 
of path-following and homotopy continuation, we show that it is no loss to assume that this 
positive chamber A is bounded. Suppose that pi{y), . . . ,pi+n{y) are degree 1 polynomials 
in y e R' that span the linear space of degree 1 polynomials, let (Pij) E ]R('+")^' be an array 
of real numbers, and consider the system of master functions 

l+n 

(1.9) n^^^^)^" = 1 j=h---,i, 

i=l 

in the polyhedron 

A := {y eM} \piiy)>0, t = l,...J + n}. 

Suppose that A is unbounded. By our assumption on pi, . . . ,Pi+n, A is strictly convex 
and therefore has a bounded face, F. Then there is a degree 1 polynomial, q{y) := ao + a ■ y 
where a G and > 0, that is strictly positive on A, and such that F is the set of points 
of A where q{y) achieves its minimum value on A. Dividing by ao, we may assume that the 
constant term of q{y) is 1. 

Consider the projective coordinate change y ^ y, where 

(1.10) yj := yj/q{y) for j = 1, . . . , / . 
Then we invite the reader to check that 

Vj = Vj/riy) for j = 1,...,/, 

where r(y) := 1 — a ■ y. Note that q{y)r{y) = 1. 

The coordinate change (11.101) manifests itself on degree 1 polynomials as follows. If p{y) 
is a degree 1 polynomial, then 



(1.11) p{y) = p{y)/r{y) and p{y) = p{y)/q{y) , 
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where p{y) := p{y) — p{0)a ■ y. Under the projective transformation (ll.lOp . the polyhedron 
A is transformed into A, where 

A := {y eR^ \ r{y) > and Pi{y) > for i = 1, ...,/ + n} . 

Proposition 1.3. Under the coordinate change fll.lOl) . the system of master functions f ll.9p 
on A is transformed into the system 

l+n 

Y[Pi{yf''' = 1 for J = 1,..., I, 

i=0 

on the polyhedron A , where Pq(jj) =q{y) and poj := — Aj. Furthermore, A is bounded. 

Proof. If y G A and y E A are related by the coordinate transformation (11.101) . then 
Pi{y) =Pi{y)/r{y)- Thus 



l+n /l+n \ l+n 

P 



i=l \i=l J i=0 

where Poijj) = r{y) and Poj = This proves the first statement. 

All that remains is to show that A is bounded. The polyhedron A is the Minkowski sum 
(1.12) A = P + R, 

where P is a bounded polytope and i? is a strictly convex polyhedral cone with vertex the 
origin. The bounded faces of A are faces of P. In particular, the bounded face F consisting 
of points where q achieves its minimum on A is also a face of P. Additionally, q achieves it 
minimum value on R at the origin. In particular, if 7^ G -R, then a-v = q{v) — q{0) > 0. 

Let 71 be an upper bound for the ratio |||/||/q'(l/) for y E P, which bounds when y 
corresponds to a point in P. Letting v range over a set of generators of the (finitely many) 
extreme rays of R shows that there is a positive constant p such that 

a ■ V > \\v\\/ p for V E R. 

Let y E A. Then y := y/r{y) is the corresponding point of A. Writing y = yp + yn where 
yp E P and yp E R, then q{y) = q{yp) + a ■ yp with both terms non-negative, and we have 



\y\ 



y 



^ \\yp\\ I Ibfill ^ llyp|| ^ Ibflll ^ 
q{y) q{y) ~ q{yp) a-yn ~ 



qiy) 

This shows that A is bounded. □ 

We remark that this coordinate transformation y \-^y is a concrete and explicit projective 
transformation transforming the unbounded polyhedron A into a bounded polyhedron A. 

1.2. Curve tracing. Let (7: M" — W^~^ be a smooth function such that 5'^"'^(0) is a smooth 
curve C. Given a point c G M" on or near C, to move along the arc of C near c is to produce 
a series of approximations to points along this arc. Predictor-corrector methods [2] are the 
standard general numerical method for this task. They proceed as follows: 
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I Prediction . Move some specified distance (the steplength) At along tlie tangent line to 
C at c. (Tlie tangent direction is the kernel of the Jacobian matrix of g at C.) That is, if 
T is the unit tangent vector to C pointing in the appropriate direction, set Cp := c + At ■ T. 

II Correction . Choose a hnear equation L{z) = vanishing at Cp and transverse to T, 
for example L{z) := T ■ {z — Cp). If At is sufficiently small, then Cp is close to a point on 
C where L = 0. Adding L to g gives a system g with Cp close to a solution, and Newton 
iterations starting at Cp for the system g approximate this nearby point on C with L = 0. 

Since C is smooth and At small, the Jacobian matrix of g has full rank in a cylinder 
around C near the point c and Cp lies in that cylinder. This procedure will produce a 
sequence of approximations Cp, Ci, . . . , to a point on C. 

III Adaptive Steplength. If correction fails to converge after a preset number (say k = 2) 
of iterations, then reduce At by some factor, say |, and return to Step I. Otherwise, if several 
consecutive predictor-corrector steps have succeeded, increase At by some factor, say 2. In 
this way, mild curves are tracked quickly and wild curves are tracked cautiously. The entire 
predictor-corrector procedure is terminated as a failure if At falls below a threshold which 
is a function of the user's patience and the available numerical precision. 

IV Resetting or Stopping. If the step was a success, decide whether to stop, based on 
stopping criteria specific to the problem being solved. If it is not time to stop, set c = Ck 
and go to Step I. If the step was a failure, leave c as is and go to Step I. 

If another branch of the curve C is near the arc we are tracing and if the steplength At 
is large enough, then Cp may fall in the basin of convergence of the other branch. The only 
general remedy for this curve- jumping is to restrict both the maximum allowed steplength At 
and the number of Newton corrections, typically to 2. Even then, curve-jumping may occur, 
only to be discovered with problem-specific checking procedures as described in Section |H 

There are many possible variations. For example, the tangent predictor could be replaced 
(after one or more initial steps) with an interpolating predictor. For more details, see [2]. 

1.3. Homotopy continuation. Predictor-corrector methods are the method of choice for 
homotopy continuation. They involve tracking a curve with a specific structure defined in 
a product of C" (the space of variables, treated as M^"") and C (a parameter space). 
For the purposes of our discussion, let 

f{z) = , with / : C" ^ C" , 

be the target polynomial system for which the isolated complex solutions are sought. Based 
upon the structure of the system /, a start system 

g{z), with ^ : C" — ^ C" , 

is chosen and the homotopy 

H{z, t) := f{z) ■ (1 - t) + 7 ■ t ■ g{z), with H : C" x C ^ , 

is formed. Here, 7^ 7 G C is a randomly chosen complex number. The start system g{z) 
is chosen so that it is easily solved and so that its structure is related to that of /. 

At t = 1, we know the solutions of H{z, 1) = 'jg{z), and we seek the solutions of H{z, 0) = 
f{z). For any given value of t, H{z,t) is a polynomial system with some number Nt of 
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isolated solutions. This number is constant, say A^, for all except finitely many t G C. For 
a general 7, restricting H{z, t) to the interval (0, 1] gives N real arcs, one for each solution 
to g{z) = when t = 1. 

Some arcs will extend to t = and their endpoints give all isolated solutions of the target 
system f{z) = 0. Some may diverge (wasting computational time), and two or more arcs 
will converge to each singular solution. Different choices of start system g{z) = will have 
different numbers of the extraneous divergent arcs. This is illustrated in Figure [H 






1 

No extraneous arcs 



1 
Singular arcs 



1 
Divergent arc 



Figure 2. Arcs in a homotopy 

Here, we track paths rather than trace curves. The predictor step is the same as in curve- 
tracing, but rather than correct along a line normal to the tangent, we freeze the value of t 
and correct back to the curve for that value of t. For more details, see [18]. 

1.4. Curve-tracing versus path-tracking. The primary difference between path-tracking 
in homotopy continuation and general curve-tracing as described in Section 11.21 is that in 
curve-tracing there are only application-specific safeguards against curve-jumping. Homo- 
topy continuation has a well-developed theory to avoid path-crossing. For almost all 7 
(those that do not satisfy a certain algebraic relation), all paths will remain distinct until 
t = 0. As a result, all paths may be tracked to some small value of t and compared to check 
for path-crossing. 

A second difference involves methods to follow ill-conditioned curves. While the curves 
to be tracked are generally non-singular, they are not necessarily well-conditioned. These 
algorithms involve repeatedly solving linear systems Ax = b, typically when A is some 
Jacobian, and the efficacy of this step depends upon the condition number k{A) of A. 
Wilkinson [20] showed that we can expect no more than P — \og{K{A)) digits of accuracy 
when solving Ax = b with P digits of precision. For homotopy continuation, adaptive 
precision techniques can overcome these potential losses in accuracy [5J. While this could 
likely be extended to general curve-tracing, the details have not yet been worked out. 

Finally, when tracing a curve with large curvature and large initial steplength At, it will 
take some time for At to become small enough so that Newton's method can succeed. There 
are many other situations in which curve-tracing will experience difficulties. However, the 
paths tracked in homotopy continuation are generally very well-behaved with gentle shifts 
in curvature, except possibly near t = 0. For that, endgames [T2| have been developed to 
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handle singularities at t = 0. Methods for handling high curvature, perhaps by adapting 
endgames to the general curve-tracing setting, are a topic for future research. 

2. The Khovanskii-Rolle algorithm 

We describe the structure of the Khovanskii-Rolle algorithm in the simple case of finding 
positive solutions in Section[2]2]and in the general case of finding real solutions in Section [2l3l 
Section H] will discuss our actual implementation of the algorithm. Some details of various 
subroutines — most notably a method for following curves beginning at singular points in 
the boundary of A — will be discussed there. 

The algorithm has two elementary routines: curve-tracing (as described in Section [1.21) . 
and polynomial system solving for a precomputation. We shall treat both elementary rou- 
tines as function calls in our description of the Khovanskii-Rolle continuation algorithm. 
The algorithm proceeds iteratively with the starting points for curve-tracing being either 
the output of the previous stage of the algorithm or solutions to certain polynomial systems 
restricted to faces of the domain polytope, A. 

We first describe the main continuation routine. This is a curve-tracing algorithm that 
uses solutions to one system of equations as input to find the solutions to another. We then 
describe our algorithm for finding all solutions to a system of master functions within the 
positive chamber, followed by the general method for all real solutions. 

2.1. Continuation step. Suppose that we have smooth functions /i, . . . , fi-i,g on a do- 
main A C with finitely many common zeroes in A. Also suppose that the zero locus of 
/i, . . . , fi-i is a smooth curve C in A. Let J be the Jacobian determinant J-det(/i, . . . , fi-i, g) 
of the functions /i, . . . , fi-i,g. 

The goal of this step is to trace C to find all common zeroes of /i, . . . , fi-i, g in A, taking 
as input the zeroes of J on the curve C and the points where C meets the boundary of A. 
This uses an extension of RoUe's theorem due to Khovanskii. We only use a very simple 
version of this theorem, which may be deduced from the classical Rolle theorem. 

Proposition 2.1 (Khovanskii-Rolle Theorem). Between any two zeroes of g along an arc 
of C there is at least one zero of J. 

Figure [3] illustrates the Khovanskii-Rolle Theorem when 1 = 2. Note that the Jacobian 
determinant J here is the exterior product df A dg. 

Khovanskii-Rolle continuation begins with a solution of J = on C and traces C in both 
directions until either it encounters a solution of (7 = on C or some other stopping criterion 
is satisfied, as described below. 

This clearly yields all solutions to (7 = on compact arcs of C. For non-compact arcs 
that have two or more solutions, the Khovanskii-Rolle Theorem guarantees that solutions 
of J = will be adjacent to solutions of g = 0. To ensure that all zeroes of (7 = on C 
are found, the algorithm also traces non-compact components of C starting from the points 
where C meets the boundary of A. Here is a more complete pseudocode description. 

Algorithm 2.2 (Continuation Algorithm). 

Suppose that /i, . . . , are different iable functions defined on a domain A C M' that 
define a smooth curve C, and g is another different iable function with finitely many simple 
zeroes on C. Let J := J-det(/i, . . . , fi-i, g) be the Jacobian determinant of these functions. 
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Figure 3. df A dg{a) < 0, df A dg{c) = 0, and df A dg{b) > 0. 

Input: Solutions S* of /i = ■ ■ ■ = fi^i = J = in the interior of A and the set T of points 
where C meets the boundary of A. 

Output: All solutions U to /i = ■ ■ ■ = fi-i = g = in the interior of A. 

For each s in S, follow C in both directions from s until one of the following occurs. 

(1) A solution u E U to = is found, 

(2.1) (2) A solution s' G S* to J = is found, or 

(3) A point t & T where C meets the boundary of A is found. 

For each t in T, follow C in the direction of the interior of A until one of the stopping 
criteria (12.11) occurs. Each point u E U is found twice. 

The computation of the sets 5* and T of solutions is presented after the full Khovanskii- 
RoUe algorithm. 

Proof of correctness. By the Khovanskii-RoUe Theorem, if we remove all points where J = 
from C to obtain a curve C°, then no two solutions of g = will lie in the same connected 
component of C°. No solution can lie in a compact component (oval) of C°, as this is also 
an oval of C, and the number of solutions = on an oval is necessarily even, so the 
Khovanskii-RoUe Theorem would imply that J = has a solution on this component. 

Thus the solutions of g = lie in different connected components of C°, and the boundary 
of each such component consists of solutions to J = on C and/or points where C meets 
the boundary of A. This shows that the algorithm finds all solutions g = on C, and that 
it will find each twice, once for each end of the component of C° on which it lies. □ 

Each unbounded component of C has two ends that meet the boundary of A. Writ- 
ing V{fi, . . . , fi-i, g) for the common zeroes in A to /i, ...,/;_!, (7, and ubc(C) for the 
unbounded components of C, we obtain the inequality 

(2.2) 2#V(/i,...,/z_i,^7) < #paths followed = 2#ubc(C) + 2#V(/i, . . . , /^^i, J) , 

as each solution is obtained twice. This gives a bound on the number of common zeroes of 
/i, . . . , fi-i, g in A and it will be used to estimate the complexity of this algorithm. 
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2.2. Khovanskii-Rolle algorithm for positive solutions. The Khovanskii-Rolle algo- 
rithm yields all solutions to a system of equations of the form 

l+n 

(2.3) ^j{y) := UpM^"-^ = 1, for j = 1, . . . , / , 
in the polyhedron 

(2.4) A := {yeR'\p,{y)>0}, 

which we assume is bounded. Here Po{y) , ■ ■ ■ , Pi+n{y) are degree 1 polynomials that are 
general in that no / + 1 of them have a common solution. Also, jSij are real numbers that 
are sufficiently general in a manner that is explained below. 

We first make some observations and definitions. In A, we may take logarithms of the 
master functions in (12. 3p to obtain the equivalent system 

l+n 

(Pj{y) := logijjiy) = ^Pi,j^ogpi{y) = forj = l,...,L 

The Khovanskii-Rolle algorithm will find solutions to these equations by iteratively applying 
Algorithm 12.21 For j = I, /—I, . . . , 2, 1 define the Jacobian determinants 

Jj := J-det(v3i, . . . , ipj] Jj+i, .. .,Ji) . 

Also, write the rational function ipj{y) = 'ipj~{y)/4'j'{y) as a quotient of two polynomials ipj' 
and ipj- Writing fj{y) for the difference ip~j{y) —ip~{y), 'ipiiy) = 1 is equivalent to fj{y) = 0. 

Algorithm 2.3 (Khovanskii-Rolle Continuation). 

Suppose that ipj and Jj for j = 1, . . . , / are as above. For each j = 0, . . . , /, let Sj be the 
solutions in A to ^ ^ 

^Pi = ■ ■ ■ = = Jj+i = ■ ■ ■ = Jl = , 
and for each j = 1, let Tj be the set of solutions to 

(2.5) fi = ■■■ = fj-i = and J^^, = . . . = J^ = Q 
in the boundary of A. 

Input: Solutions 5*0 of Ji = ■ ■ ■ = J; = in A, and sets Ti, . . . , T;. 
Output: Solution sets Si, ... , Si. 

For each j = 1, . . . , /, apply the Continuation Algorithm 12.21 with 

(/l5 • • • 5 fl-l) = (</5l5 • • • 5 5 Jj + l, ■ ■ ■ ,Jl) 

and g = ipj. The inputs for this are the sets Sj-i and Tj, and the output is the set Sj. 

The last set computed. Si, is the set of solutions to the system of master functions (12. 3 p 
in A. We describe our implementation for positive solutions when / = 2 in Section |4] and 
illustrate the algorithm in Section [31 

Proof of correctness. Note that for each j = 1, . . . ,1, the jth step finds the solution set Sj, 
and therefore preforms as claimed by the correctness of the Continuation Algorithm 12. 2[ 
The correctness for the Khovanskii-Rolle algorithm follows by induction on j. □ 
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The Khovanskii-Rolle algorithm requires the precomputation of the sets 5*0 and Ti, . . . , T/. 
The feasibility of this task follows from Lemma 3.4 in which we state below. By our 
assumption on the generality of the polynomials pi, a face of A of codimension k is the 
intersection of A with k of the hyperplanes Pi{y) = 0. 

Proposition 2.4. Let pi,i = 0, . . . , l+n, ipj, Jj,j = 1, . . . ,1, and A be as above. 

(1) The solutions Tj to f l2.5l) m the boundary of A are the solutions to 

Jj+i = ■ ■ ■ = Jl = Q 
in the codimension- j faces of A. 

l+n 

(2) Jj ■ Y\_Pi{y)'^ ^ polynomial of degree 2'"%. 

1=0 

By Proposition 12.4( 1). solving the system fl2.5l) in the boundary of A is equivalent to 
solving (at most) ( ^"^"^) systems of the form 

(2.6) Pi, = • • • = pi^ = Jj+i = ■ ■ ■ = Jl = 

and then discarding the solutions y for which Pi{y) < for some i. Replacing each Jacobian 

Jj by the polynomial Jj '■= Jj ■ Y[\=o PiivY' ^ ^^'^ using Pi^{y) = ■ ■ ■ = pi- = Q io eliminate 
j variables, we see that fl2.6p is a polynomial system with Bezout number 

n-2n-4n--- = 2^''^')n^-^ 

'ml — j variables. Thus the number of solutions Tj to the system (12. 5p in the boundary of 
A is at most 2C2')n'-J ('+^+^) . 

Likewise 5*0 consists of solutions in A to the system fl2.6p when j = and so has at most 

the Bezout number 2^^n^ solutions. Thus the inputs to the Khovanskii-Rolle Algorithm are 
solutions to polynomial systems in / or fewer variables. 

Theorem 2.5. The Khovanskii-Rolle Continuation Algorithm finds all solutions to the 
system (12.30 in the bounded polyhedron A (12. 4p . when the degree 1 polynomials Pi{y) and 
exponents f3ij are general as described. It accomplishes this by solving auxiliary polynomial 
systems (12. 6p and following implicit curves. For each j = 0, . . . , l—l, it will solve at most 
^i+n+i^ pglyjiQjjiidl systems in l—j variables, each having Bezout number 2^ 2^)77,'-^. Jn all, 
it will trace at most 

(2.7) /2(2)+V + ^(/+l-j)2'-V-^CY') < /^^2(2)n' 
implicit curves in A. 

Proof. The first statement is a restatement of the correctness of the Khovanskii-Rolle Con- 
tinuation Algorithm. For the second statement, we enumerate the number of paths, follow- 
ing the discussion after Proposition 12. 4[ Let Sj^tj be the number of points in Sj and T,, 
respectively, and let be the number of paths followed in the jth step of the Khovanskii- 
Rolle Continuation Algorithm. 
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By ( 12.21) . rj = tj + 2sj_i and Sj < \tj + Sj_i. So rj<tj + --- + ti + 2so, and 



ri + ---+n < 2so + . 




(2.8) So < 2(2)^' and tj < 2C^')n'-^ f+"+^) , 



we obtain the estimate on the left of (12.71) . Since l+l—j < I, we bound it by 




which is bounded by / ^^2^'^n\ by Lemma 3.5 in [7]. 



□ 



Remark 2.6. The bound (12. 7p on the number of paths to be followed is not sharp. First, 
not every system of the linear polynomials 



defines a face of A. Even when this defines a face F of A, only the solutions to (12.61) that 
lie in F contribute to T,-, and hence to the number of paths followed. Thus any slack in 
the estimates (12.80 reduces the number of paths to be followed. Since these estimates lead 
to the fewnomial bound for s/, we see that the Khovanskii-RoUe Continuation Algorithm 
naturally takes advantage of any lack of sharpness in the fewnomial bound. 

This may further be improved if A has m < I + n + 1 facets for then the binomial 
coefficients in (12.71) become (^). 

2.3. Khovanskii-Rolle algorithm for all real solutions. Section [2l2] describes how to 
find solutions to a system of master functions (12.31) in the polyhedron A (12. 4p . which is 
assumed bounded. Through Gale duality and the coordinate transformation (ILIOI) . this 
gives a method to find all positive real solutions to a system of polynomial equations (11.61) . 

To find all real solutions to a system of polynomial equations or of master functions, one 
could simply repeat this process for every chamber in the complement of the hyperplanes 
Pi{y) = for i = 1, . . . , l+n. This is however inefficient as our method (homotopy continu- 
ation) for computing the sets 5*0 and Tj for i = 1, ... ,1 of starting points for one chamber 
gives the starting points for all chambers. Besides careful bookkeeping, this requires some 
projective coordinate transformations so that the tracking occurs in bounded chambers. 

Each point in Tj is incident upon 2^ chambers, and thus is the starting point for 2^ arcs 
to be followed in Algorithm 12.31 Surprisingly, this has a mild effect on the complexity, 
requiring only that we replace the in (12.71) by e^. This observation, which was made 
while developing this Khovanskii-Rolle algorithm, was the genesis of the bound in [5]. 



Phiy) 



Pijiy) = 0, 



3. Examples 



We first illustrate the Khovanskii-Rolle algorithm and our implementation on the master 
function system of Example II. 1[ Write the system (II. 3p of master functions in logarithmic 
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form, (pi(y) = v^2(y) = 0, for y E A, which is the heptagon depicted in Figure [H Here, we 
have / = 2 and n = 5 with 7 hnear polynomials. 

ifiiy) = - \og{y) + log(2x+2y-l) - 2 log(4x) + log(i±f^) 

-21og(5^) +21og(^^^) -log(3-2x+i/). 

ip2{y) = log(?/) + 61og(2x+2?/-l) -31og(4x) + 61og(^i±f^) 

-21og(^^) + 71og(«^^) + log(3-2x+y) . 

The polynomial forms J2, Ji of the Jacobians are (omitting the middle 60 terms from Ji), 

J2 = -168x^ - mex'^y + mx^^ - 536x^y^ - 1096x?/^ + 456?/^ + 1666a;^ + 2826x^?/ 
+3098x^y^ + 690Axy^ - 1638?/^ - 3485a;3 - 3721a;2y - 153182;?/2 _ 1836^^ 
+1854x2 + 8AA2xy + 9486?/^ - 192x - 6540?/ + 720 . 

Ji = 10080x^° - 168192x^1/ - 611328xV - ■■■ + 27648x + 2825280?/ . 

Remark 3.1. Instead of the polynomial form J2 of the Jacobian of cpi and ip2, we could 
use the Jacobian J{f,g) of / and g (11 ■4p . This however has degree 25 with 347 terms and 
coefficients of the order of 10^''. The Jacobian of this and / has degree 29 and 459 terms. 
This control on the degree of the Jacobians is the reason that we use the logarithms of the 
master function in the formulation of the Khovanskii-RoUe Algorithm. 

We now describe the Khovanskii-RoUe algorithm on this example. 
Precomputation. We first find all solutions Sq to Ji = J2 = in the heptagon A, and all 
solutions J2 = in the boundary of the heptagon. Below are the curves J2 = and Ji = 
and the heptagon. The curve J2 = consists of the four arcs indicated. The remaining 
curves in this picture belong to Ji =0. On the right is an expanded view in a neighborhood 
of the lower right vertex, (|,0). 




A numerical computation finds 50 common solutions to Ji = J2 = with 26 real. Only six 
solutions lie in the interior of the heptagon with one on the boundary at the vertex (3/2, 0). 
There are 31 points where the curve J2 = meets the lines supporting the boundary of 
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the heptagon, but only eight he in the boundary of the hexagon. This may be seen in the 
pictures above. 

First continuation step. Beginning at each of the six points in A where Ji = J2 = 0, 
the algorithm traces the curve in both directions, looking for a solution to '■Pi = J2 = 0. 
Beginning at each of the eight points where the curve J2 = meets the boundary, it follows 
the curve into the interior of the heptagon, looking for a solution to v^i = </2 = 0. In tracing 
each arc, it either finds a solution, a boundary point, or another point where Ji = J2 = 0. 
We may see that in the picture below, which shows the curves ipi = Q and J2 = 0, as well 
as the points on the curve J2 = where Ji also vanishes. 




In this step, 2 ■ 6 + 8 = 20 arcs are traced. The three solutions of J2 = v^i = will each be 
found twice, and 14 of the tracings will terminate with a boundary point or a point where 
Jl = 0. 

Second Continuation step. This step begins at each of the three points where v^i = J2 = 
that were found in the last step, as well as at each of the six points where = meets 
the boundary of the heptagon (necessarily in some vertices). 




Curve tracing, as described in more detail in the next section, can be carried out even 
in the presence of singularities, as in the case of the curves initiating at vertices. In this 
case, this final round of curve-tracing revealed all six solutions within the heptagon, as 
anticipated. Furthermore, each solution was discovered twice, again, as anticipated. 
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By Theorem I2.5[ the bound on the number of paths followed (using 7 = / + n in place of 
/ + n + 1 in the binomials as in Remark 12.61) is 

2 ■ 2 ■ 2(2) + 2 ■ 2(2) -5^ ■ (l) + 1 ■ 2(2) ■ 5° ■ Q = 298 . 
By Theorem 3.10 in [7J, the fewnomial bound in this case is 

2.52^^ (5+l)(5+3) j ^ 

In contrast, we only traced 20 + 12 = 32 curves to find the six solutions in the heptagon. 
The reason for this discrepancy is that this bound is pessimistic and the Khovanskii-RoUe 
Continuation Algorithm exploits any slack in it. 

Timings and comparison to existing software. The system of Laurent polynomi- 
als (11.11) was converted into a system of polynomials by clearing denominators. It was then 
run through PHCpack, Bertini, and the proof-of-concept implementation described in the 
next section. As described in Section [H this system has 102 regular solutions, 10 of which 
are real. All runs of this section were performed on a 2.83 GHz running CentOS with Maple 
13, Bertini 1.1.1, and PHCpack v2.3.48. 

In blackbox mode, PHCpack used polyhedral methods and found 102 regular solutions in 
around 2.5 seconds, though it only classified eight of them as real. Using all default settings 
and a 5-homogeneous start system, Bertini found all 102 regular solutions and identified 
the 10 that are real. However, because of the use of adaptive precision, this took around 
23 seconds. Using fixed low precision, this time dropped to around 9 seconds while still 
identifying the solutions correctly. 

The implementation described in the next section found all positive solutions of the Gale 
dual system of master functions 11.41 in around 15 seconds using safe settings. Almost all of 
this time was spent in computing the sets 5*0 and Ti in Bertini, using adaptive precision for 
security. By changing to fixed low precision for the Bertini portions of the computation, the 
timing for the entire run fell to around 1.4 seconds - the shortest time of all runs described 
here. Though such efficiency is welcome, security is more valuable. It is expected that more 
sophisticated software than that described in the next section will be more efficient. 
A more extreme example. A polynomial system with high degree, many complex solu- 
tions, and few real solutions further illustrates the value of Khovanskii-RoUe continuation. 
For example, consider the system of Laurent polynomials 



10500 - tu^^^ - 3500t-^u^^^v^w^ = 

. 10500 - t - 3500riM69Vw5 = 

14000 - 2t + _ 3500^, = 

14000 + 2t - _ 35oo^y = 0. 

By solving a set of master functions Gale dual to (13. ip for 0, we obtain 



3500^23,8^4(3 _^)45_ (^3 _^)33(4_2a; + ^)60(2a;_^+ 1)60 ^ g 
3500^23,27(^3 _3,)8(3_^)4_ ^15(4 _ + ^)60j2a;_^+ 1)60 ^ Q_ 
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System (13. ip has 7663 complex solutions but only six positive real solutions. PHC com- 
putes these solutions in 39 minutes, 38 seconds while the implementation of the next section 
takes only 23 seconds to find them, using safe settings. 

4. Implementation 

We have implemented the Khovanskii-Rolle Algorithm to find all solutions to a system of 
master functions in a bounded polyhedron A, but only when / = 2. This proof-of-concept 
implementation relies on external calls to the Bertini software package [4j as a polynomial 
system solver. The implementation is in a Maple script that is publicly available at We 
plan to implement a general version of the Khovanskii-Rolle algorithm in the future. 

We describe some aspects of this implementation, including the precomputation to find 
the solution sets Sq and Ti and the curve-tracing method from these points and from 5*1, all 
of which are smooth points on the traced curves. We also discuss tracing curves from the 
vertices T2, which is non-trivial as these curves are typically singular at the vertices. Lastly, 
we discuss procedures for checking the output. 

4.1. Polynomial system solving. The precomputation of Sq and Ti for Algorithm 12.31 
involves finding the solutions of systems of Jacobians (12.61) within A and on its boundary. 
For this, we use the system Bertini For each system to be solved, the maple script 
creates a Bertini input file, calls Bertini, and collects the real solutions from an output file. 
Bertini, as with all homotopy methods, finds all complex solutions. However, as explained 
previously, this overhead may be much less than would be encountered in the direct use of 
homotopy methods to solve the original system. 

4.2. Curve-tracing from smooth points. The points of Sq and Ti from which we trace 
curves in the first step of the algorithm, as well as the points 5*2 used in the second step, are 
smooth points of the curves J2 = and (fi = 0, respectively. This curve-tracing proceeds as 
in Section 11.21 Indeed, suppose we are tracing a curve C in the polytope A, starting from 
points where C meets the boundary of A and from interior starting points where a Jacobian 
determinant J vanishes, and we seek points of C where some function g vanishes. Then, as 
described in the Continuation Algorithm 12. 2[ there are three basic stopping criteria: 

(1) The tracer passes a point where g = 0, which is a point we are looking for. 

(2) The tracer passes a point where J = 0, which is another starting point. 

(3) The tracer leaves the polytope. 

The computation of the tangent and normal lines is straightforward, as is the linear 
algebra required for curve-tracing. Each predictor-corrector step gives a point p near C. 
We compute the values of g, J, and the degree 1 polynomials Pi defining A at p. A sign 
change in any indicates a stopping criteria has been met. When this occurs, bisection in 
the step length At from the previous approximation is used to refine the point of interest. 

It may seem that these (sign-based) criteria will fail when there are clustered (or multiple) 
solutions to (7 = on the curve between the current and previous approximation. However, 
the Khovanskii-Rolle Theorem implies there will be zeroes of the Jacobian determinant J 
interspersed between these solutions (or coinciding with the multiple solutions). Further- 
more, the number of solutions to g = and to J = will have different parities, so that 
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one of the two functions g and J will change sign, guaranteeing that one of the stopping 
criteria will be triggered for such a curve segment. 

Curve-tracing involves a trade-off between security and speed. Security is enhanced with 
a small steplength At and allowing only one or two Newton correction steps. However, these 
settings contribute to slow curve-tracing. The examples in Section [3] were computed with 
secure settings. In Section 14.41 we give methods to detect some errors in this curve-tracing 
which can allow less secure settings. 

4.3. Curve-tracing from vertices. The curves to be traced are typically singular at the 
vertices of A. While traditional curve-tracing fails at such points, we employ a simple 
alternative that takes advantage of the local structure of the curves. 

At a vertex v of A, we make a linear change of coordinates so that the two incident edges 
are given by the coordinate polynomials yi and y2- Then (pi{y) = may be expressed as 

n+2 

llMyf-' = 1. 

j=0 

Setting Pn+i = yi and Pn+2 = 2/2, we may solve for ?/2 to obtain 



2/2 = Z/l 

i=0 



In the neighborhood of v this is approximated by the monomial curve 

n 
=0 



where we have evaluated the terms Pi{y) for i < n at the vertex v, i.e., the product is just 
a constant a. We may write this expression as 

(4.1) y2 = a-yf ; 

note that a > 0. 

We do not begin tracing the curve fi{y) = at t>, but instead begin from a point Cp on 
the monomial curve (14.11) in A near v. If the first predictor-corrector step from Cp succeeds 
to approximate (pi{y) = 0, then we trace this curve as described before. If this predictor- 
corrector step fails from Cp, then we simply choose a point on the monomial curve a bit 
further from v and try again. In our experience, this special form of monomial tracking 
quickly gives way to usual curve-tracing on (pi{y) = 0. 

4.4. Procedures for checking the output. General curve-tracing is not foolproof. If 
the curves to be traced are very close together, then curve-jumping may occur, leading to 
missed solutions. However, for curve-tracing in the Khovanskii-RoUe algorithm, there is a 
simple way to check for such errors. 

As described after Algorithm 12. 2^ each point in the sets St for i = should be 

discovered twice as the end of an arc that is tracked. Thus, if a solution is not found twice, 
an error occurred in the curve-tracing. 

This check will not capture all errors. The development of further verification and certi- 
fication procedures is a goal for future research. This lack of checks is not uncommon for 
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new algorithms, including the methods introduced in the early development of numerical 
algebraic geometry 



5. Conclusions 

Numerical homotopy continuation is a robust and efficient algorithm for finding all com- 
plex solutions to a system of polynomial equations. Real solutions are obtained by selecting 
those solutions with small imaginary parts. While often practical, this is wasteful and does 
not exploit the real algebraic nature of the problem. 

We presented a numerical continuation algorithm to find all real or positive solutions to 
a system of polynomials, along with details of our implementation and two examples. This 
new Khovanskii-RoUe continuation algorithm is efficient in that the number of paths to be 
followed depends upon the corresponding fewnomial bounds for the numbers of real solutions 
and not on the number of complex solutions. This is a significant difference between our 
new algorithm and all other known methods for solving polynomial systems. 

This algorithm does not directly solve the given polynomial system but rather an equiva- 
lent (Gale dual) system consisting of master functions in the complement of an arrangement 
of hyperplanes. This appears to be the ffist curve-tracing algorithm that finds only real so- 
lutions to a system of equations. 

This paper is the ffist step in this new line of research. There are clear generalizations 
to higher dimensions. We plan further research into techniques for tracing the curves that 
begin at singularities. In contrast to path-following in homotopy continuation, the security 
of curve-tracing in this algorithm relies on heuristics. Another research direction is to 
enhance the security and efficacy of curve-tracing. 

When our algorithm has been implemented in more than two variables (Z > 2), we plan 
to use it to study real solutions to systems of polynomial equations. 
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